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ABSTRACT: The direct detection of gravitational wave by Laser Interferometer Gravitational- 
Wave Observatory indicates the coming of the era of gravitational-wave astronomy and 
eravitational-wave cosmology. It is expected that more and more gravitational-wave events 
will be detected by currently existing and planned gravitational-wave detectors. The grav- 
itational waves open a new window to explore the Universe and various mysteries will 
be disclosed through the gravitational-wave detection, combined with other cosmological 
probes. The gravitational-wave physics is not only related to gravitation theory, but also 
is closely tied to fundamental physics, cosmology and astrophysics. In this review article, 
three kinds of sources of gravitational waves and relevant physics will be discussed, namely 
gravitational waves produced during the inflation and preheating phases of the Universe, 
the gravitational waves produced during the first-order phase transition as the Universe 
cools down and the gravitational waves from the three phases: inspiral, merger and ring- 
down of a compact binary system, respectively. We will also discuss the gravitational waves 
as a standard siren to explore the evolution of the Universe. 
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1 Introduction 


On 11 February 2016, the LIGO Scientific Collaboration and the Virgo Collaboration [1] 
announced that on 14 September 2015 at 09:50:45 UTC the two detectors of the Laser In- 
terferometer Gravitational-Wave Observatory (LIGO) simultaneously observed a transient 
gravitational-wave (GW) signal. The GW event is named GW150914. The frequency of 
the signal increases from 35 to 250 Hz with a peak GW strain 1.0 x 107?! at the 150 Hz. 
The GW signal is consistent with the one predicated by general relativity for the inspiral 
and merger of a pair of black hole and the ringdown of the resulting single black hole. 
The source is located at a luminosity distance of 4104460 Mpc corresponding to a redshift 
Z= 0.094993. The initial black hole masses are 3615 Mo and 2914 Mo and the resulting 
black hole mass is 62+4 Mo. The mass difference 3.0+85 Mo is radiated in the form of GWs. 
This is the first direct detection of GW and the first observation of a binary black hole 
merger. On 15 June 2016, the second GW event, GW151226, was announced by the same 
team [|2]. This time, observed signal lasts approximate 1 s, the frequency increases from 
35 to 450 Hz within 55 cycles, and the peak gravitational strain reaches ZAER x 107-22 
at 450 Hz. The source is also the merger of two black holes with masses 14.2+33 M. and 


7.5423 Mo, respectively, the final black hole has mass 20.8797 Mo. This event happened 
with a luminosity distance 4401180 Mpc corresponding to a redshift 0.09+9:03, 

The GW was predicted by Albert Einstein in 1916 [3, 4], 1 year later after he finally 
formulated his theory on gravitation, genera relativity. But the physical reality of the GW 
solution of the Einstein field equations was not showed until the Chapel Hill conference in 
1957 [5]. In [6, 7], it has been shown that GW carries energy and when passing through 
the spacetime in a form of a sandwich, it affects test particles. More than one century has 
passed since the Einstein’s proposal of general relativity, although it passed various precise 
tests, some alternatives still survive, for example, scalar-tensor gravity theory, f (R) gravity, 
modified gravity with higher curvature terms, etc. On the other hand, based on general 
relativity, one has now a standard model (SM) of cosmology, ACDM model, which is quite 
well consistent with various astronomical observations made so far. Now we understand 
well that GW exists not only in general relativity, but also in other relativistic covariance 
gravity theories. Due to the limited space, this review is confined for the GWs in general 
relativity. 

The sources of GWs could be classified into two categories roughly. One is called cos- 
mological origin, the other is relativistic astrophysical origin. In the cosmological case, GWs 
can be produced in the early stages of the Universe, for example, during the inflation and 
reheating epochs. Such GWs are called primordial GWs, and they will leave unique imprint 
on the cosmic microwave background (CMB), the so-called B-mode. On the other hand, 
during the evolution of the Universe, it is expected that various phase transitions (PTs) 
had happened as the temperature of the Universe decreases, for example, the symmetric 
breaking of the grand unification theory, electroweak (EW) PT, quantum chromodynamics 
(QCD) PT, etc. During those PTs, GWs are also expected to be produced with different 
features. Also the interactions of topological defects such as cosmic string, domain wall, 
etc, produced during PTs, will create GWs. Therefore, the detection of GWs due to those 
cosmological origins can reveal physics associated with the evolution of the Universe. In the 
astrophysics side, GWs can be produced in various processes, for example, rotation of non- 
symmetric neutron star, explosion of supernovae, inspiral, merger and ringdown of some 
compact binaries including white dwarf, neutron star and/or black hole. In particular, the 
compact binary systems are main sources for GW detection such as LIGO, Virgo, Kagra, 
Einstein Telescope, etc., ground-based GW experiments, and Laser Interferometer Space 
Antenna (LISA), Deci-Hertz Interferometer Gravitational wave Observatory (DECIGO), 
Big Bang Observer (BBO), Taiji, Tianqin, etc., space-based GW experiments. 

Therefore, the GW physics is closely related to fundamental physics, cosmology and 
astrophysics. The direct detection of GWs in |1, 2] indicates the coming of the era of GW 
astronomy and GW cosmology. The detection of GWs opens a new window to explore the 
Universe. Combining the electromagnetic (EM) radiation, neutrino, cosmic ray and GWs, it 
could be expected that one is able to reveal various currently existing mysteries concerning 
the early evolution of the Universe, property of dark matter, the nature of dark energy, etc. 
In this brief review, we are going to summarize some important aspects relevant to the GW 
physics. 

The outline of this review is as follows. Section 2 is going to introduce the GWs 


produced in the primordial Universe and its detection through CMB. In particular, we 
emphasize the properties of GWs created during inflation and preheating processes and 
current situation of detection of the primordial GWs. In section 3 we will discuss the GWs 
produced during cosmic PT. There we will introduce the bubble nucleation, bubble expan- 
sion and bubble percolation. During the strong first-order PT, bubble collision, turbulent 
magnetohydrodynamics (MHD) and sound wave are all the sources to produce GWs. Sec- 
tion 4 will mainly be devoted to discussions on the GWs from the dynamics of compact 
binary systems. There we will introduce three main methods to solve the binary system: 
the post-Newtonian (PN) approximation, numerical relativity and the black hole perturba- 
tion, corresponding three phases: inspiral, merger and ringdown of two black holes system, 
respectively. In that section, we will also discuss the possibility of GWs produced by binary 
systems as a cosmological probe. With this new probe, it is expected to have a strong 
constraint on cosmological parameters, combining with other cosmological probes. 


2 Gravitational Waves From Primordial Universe 


In order to solve some problems in big-bang cosmology such as the horizon and flatness 
problems, inflationary scenario was introduced [8-10], in which a period of accelerated ex- 
pansion of the Universe happened at early times. Inflation not only predicts the primordial 
scalar perturbations, which provide a natural way to generating the anisotropies of the 
CMB radiation and the initial tiny seeds of the large-scale structure observed today in the 
Universe, but also generates a stochastic background of the primordial GWs. Although such 
a stochastic background of the primordial GWs has not been observed yet, its detection 
would open a new window to understanding the physics of the early Universe and thus the 
origin and evolution of the Universe. In this section, we shall firstly review the properties of 
the primordial GWs produced during inflation and preheating (see Fig. 1), and then discuss 
observational implications. 

GWs are described by a transverse-traceless gauge-invariant tensor perturbation, hij, 
in a Friedman-Robertson-Walker (FRW) metric, 


ds? = a? (T) [—dr? + (Sij + hi; )dx'da | , (2.1) 


where 7 is the conformal time and a is the scale factor, and hj; satisfies Ohi; = 0 and 
gti hij = 0. To first order in hij, the perturbed Einstein equation reads 


2 
hi, + 2Hhij — V7 hij = wets (2.2) 
P- 


where the prime denotes the derivative with respect to T, H = a’ /a is the Hubble parameter 
in T, Mp = (8G)? is the reduced Planck mass, and the source term; mt is the 
transverse-traceless projection of the anisotropic stress tensor Tj;. Since hij is symmetric, 
transverse and trace-free, tensor modes are left with two physical degrees of freedom, which 
are expanded in the Fourier space as 
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Figure 1. The schematic illustration of the inflaton potential. In the slow-roll inflationary scenario, 
an accelerated expansion of the Universe occurs when the inflaton rolls slowly along its potential. 
After inflation ends, the inflaton oscillates around the minimum of its potential, whose energy is 
converted into radiations. 
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where e are the polarization tensors with two polarization states (+, x) of GWs. We 


have the power spectrum of GWs as 
Pr(rT,k) = I ht? + |h]? 2.4 
(rk) = Eat + P, (2.4) 
and the energy spectrum of GWs as Qgw(7T, k) = dpgw/dln k/pe, where 
1 B 
M2 (Hh) (2.5) 


Pew(T, k) = Ag? 


is the energy density and pe = 3H MÀ is the critical density of the Universe. 


2.1 During Inflation 


In the standard single-field slow-roll inflationary scenario, at first order in perturbation 
theory Eq. (2.2) reduces to a free wave equation. In this case, the wave equation can 
analytically be solved in the slow-roll approximation. The Bunch-Davies vacuum condition 
in the asymptotic past is imposed because the modes lie well inside the Hubble radius. Of 
course, if a general vacuum condition is imposed, there are additional features in the power 
spectrum. During inflation, quantum fluctuations are amplified and stretched, and then 
nearly frozen on super-Hubble scales. The single-field slow-roll inflation predicts a slightly 


Pr(k) = ay (2.6) 


red-tilted spectrum: 


and a consistency relation nr = —r/8 between the tensor spectral index ny and the tensor- 
to-scalar ratio r. Here, np and r are evaluated at the epoch when a given perturbation mode 
leaves the Hubble horizon, i.e. k = aH. We usually introduce the tensor-to-scalar ratio 
r = Ar/As at a pivot scale kx, i.e. the amplitude of tensor perturbations Ar with respect 
to that of scalar perturbations Ag. The red-tilted spectrum means that the amplitude of 
tensor perturbations becomes small on small scales due to the fact that large-scale modes 
are earlier stretched across the Hubble horizon than small-scale modes as the energy density 
slowly decreases during inflation. From the wave equation we see that the evolution of hij 
depends explicitly on a(r) and implicitly on inflation potential through FRW equations. 
Hence the power spectrum of GWs encodes useful information on the evolution of the scale 
factor. Moreover, the tensor-to-scalar ratio is related to the energy scale of inflation by 
V= 3r’ My Asr/2 = (1.88 x 10'6GeV)4r/0.10. Here we have adopted the estimated value 
of the amplitude of scalar perturbations from the Planck 2015 data [11]. Different models 
predict different values of the tensor-to-scalar ratio. For example, the simplest chaotic 
inflation with a quartic potential [10] predicts a large value of r ~ 0.26 while the R? 
inflation [12] predicts a small value of r ~ 0.0033 to lowest order in slow-roll parameters if 
the number of e-folds N, = 60 is assumed. With the help of the scalar spectral index, the 
estimated value of r is robust to discriminate slow-roll inflationary models. In summary, 


the measurement of the power spectrum of GWs helps us to 
e test the vacuum initial condition, 
e detect the evolution of the scale factor, 
e determine the energy scale of inflation, 
e discriminate inflationary models. 


If the primordial GWs are detected, the next important question to answer is what 
is the shape of the power spectrum of GWs and whether there are additional features in 
the power spectrum. In the slow-roll inflationary scenario, the shape of the tensor power 
spectrum is characterized by the tensor spectral index nr since the running of the spectral 
index is negligible to the lowest order in slow-roll parameters. More general shapes beyond 
slow-roll may be reconstructed by using a binning method of a cubic spline interpolation 
in a logarithmic wavenumber space [13-15]. Checking the consistency relation between 
nr and r provides a powerful test of the single-field slow-roll inflationary scenario. The 
violation of the consistency relation could in principle come from the following two aspects. 
The first is that the second and third terms on the left hand side of Eq. (2.2) are modified 
in general inflationary models, such as the k-inflation [16], Gauss-Bonnet inflation [17— 
19] and generalized G-inflation [20]. In the model of k-inflation, the consistency relation 
becomes np = —r/8cg, where cg is the sound speed of scalar perturbations. In the Gauss- 
Bonnet inflationary model, the consistency relation is broken due to nr = —r/8 — 61, 
where 6; is determined by the Gauss-Bonnet coupling term. In the model of generalized 
G-inflation, the tensor spectral index not only depends on the evolution of the scale factor 
but also the higher-order derivative terms, which admits a blue-tilted power spectrum of 


GWs. The second is that there exists a non-negligible source term on the right hand side of 
Eq. (2.2) during inflation. In this case, the wave equation is solved by the Green’s function 
method. Possible sources of generating GWs include first-order scalar perturbations [21], 
perturbations of the extra field such as the curvaton |22] and spectator field |23, 24], and 


particle production during inflation [25]. 


2.2 During Preheating 


In the inflationary scenario, at the end of inflation, the inflaton field begins to oscillate 
around the minimum of its potential. Such coherent oscillations produce elementary parti- 
cles and eventually reheat the Universe. This process is called reheating [26]. The coupling 
between the inflaton field and other fields is necessarily tiny ensuring that reheating proceeds 
slowly. Preheating provides a more rapidly efficient mechanism for extracting energy from 
the inflation field by parametric resonance [27]. Such a process is so rapid that the produced 
particles are not in thermal equilibrium. The preheating leads to large and time-dependent 
inhomogeneities of the stress tensor that source a stochastic background of GWs [28]. Un- 
like GWs produced during inflation, they are generated and remain in the Hubble horizon 
until now. Clearly, their wavelengths are smaller than the Hubble radius at the time of 
GW production. Therefore, the peak frequency of this type of stochastic GWs is typically 
of order more than 10° Hz. Detecting such high-frequency GWs is particularly challenging. 
For example, for the ¢* and ¢? chaotic inflationary models, lattice simulations [29] show 
that preheating can lead to GWs with frequencies of around 10° ~ 108 Hz and peak power 
of Qgwh? ~ 107? ~ 1071H at present [30]. For hybrid inflation, GWs cover a larger range 
of frequencies. The peak wavelength depends essentially on the coupling constant [31]. 
See [32] for the most recent review on the GWs from the inflationary era and preheating 
era. 


2.3 Observational Implications 


Observational constraints on GWs produced during inflation mainly come from the B-mode 
polarization of the CMB anisotropies. Such GWs generate a quadrupolar anisotropy of 
momentum m = 2 in the intensity field of photons at the epoch of recombination while scalar 
perturbations generate only a quadrupolar anisotropy of momentum m = 0. Importantly, 
only the quadrupolar anisotropy of momentum m = 2 causes the B-mode polarization. 
Therefore, the measurement of the B-mode polarization allows us to probe GWs produced 
during inflation. The Planck 2015 data give the 95% confidence level upper bound for 
the tensor-to-scalar ratio ro.oo2 < 0.10 at the pivot scale ką = 0.002 Mpc™! [11], which is 
improved to ro.o02 < 0.08 by adding the BKP cross-correlation likelihood [33]. With the help 
of the scalar spectral index ng, slow-roll inflationary models are discriminated in the r—ng 
plane, as shown in Fig. 2. For example, the inflationary model with a quartic potential [10] 
is strongly disfavored by recent CMB data. Future ground-based and space-based CMB 
experiments will provide high precision measurements of the B-mode polarization. Actually, 


constraints on GWs produced during inflation from measurements of the CMB anisotropies 


are limited to a narrow scale range of 1074 ~ 1071 Mpc~!, which corresponds to the change 
of e-folds number AN ~ 8. This implies that it is impossible to detect the global shape 
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Figure 2. Marginalized joint 68% and 95% confidence level regions for ns and ro.oo2 from the 
Planck 2015 data, compared to the theoretical predictions of some inflationary models. This figure 
is quoted from [11]. 


of the tensor spectrum by using only CMB data. Laser interferometer gravitational-wave 
experiments provide the possibility of measuring the stochastic background at small scales. 
During inflation, quantum fluctuations of h;; are amplified and stretched across the Hubble 
horizon, and then nearly frozen on super-Hubble scales. After reentering the Hubble horizon 
at the epoch of radiation or matter dominance, tensor perturbations begin to oscillate with 
the amplitude damped by a factor a~!. Unfortunately, for the nearly scale-invariant tensor 
spectrum predicted by the single-field slow-roll inflationary models, the energy spectrum in 
the frequency range of 1071? ~ 10° Hz at present becomes too weak to be directly detected 
by pulsar timing array (PTA) experiments and laser interferometer experiments [34]. If it 
is a blue-tilted tensor spectrum predicted by inflationary models beyond slow-roll or caused 
by the source term, the direct detection might be possible in the future [35]. 

For GWs produced during preheating, since their wavelengths are smaller than the 
Hubble radius at the time of GWs production, the corresponding peak frequencies are 
typically of order more than 10? Hz. It is particularly challenging to detect such high- 
frequency GWs by ground-based laser interferometer experiments. 


3 The Gravitational Waves From Phase Transitions 


It is believed that our observable Universe has experienced several PTs (PTs), during which 
an unknown high-degree symmetry is subsequently relaxed down to the broken EW symme- 
try that is well described by the SM of particle physics. If the PT from a high-temperature 
symmetric phase to a low-temperature symmetry-broken phase is of first order, the true 
vacuum bubbles would be nucleated within the false vacuum, leading to the expanding, 
colliding and merging bubbles that generate a stochastic background of GWs (see [36] for 


recent review and [37| for the discussions of GWs from cosmic strings and domain walls, 
which will not be discussed in this review). The primary motivations to study the GWs 
from PTs are two-folds: First, the EWPT of SM is cross-over according to the current mea- 
surements of Higgs mass. Therefore, any detections of GWs from PTs would necessarily 
provide us a unique probe beyond the SM, which cannot be directly probed by the particle 
colliders in a foreseeable future; second, the stochastic backgrounds of GWs also consist of 
the GWs from the inflation era and the reheating era and the other cosmological defeats 
from PTs like cosmic strings and domain walls. Therefore, it would help us to extract the 
GWs of astrophysical sources from the stochastic backgrounds of GWs. However, the topics 
of detections are not discussed in this review, which should merit another paper to discuss 
how to detect the stochastic background of GWs and how to distinguish the signals of PTs 
from those signals of reheating or other cosmological defeats. 


3.1 Bubble Nucleation 


In the seminal papers [38, 39], the first semiclassical description of vacuum decay in flat 
spacetime was developed for a single scalar field without derivative interactions. The vac- 
uum decay is implemented through barrier penetration from the unstable false vacuum to 
the stable true vacuum, of which the field configuration is captured in the so-called bounce 
equation. The bounce equation describes the true vacuum bubbles nucleated during barrier 
penetration in the surrounding false vacuum with the probability per unit time per unit 
volume being of form IT = Aexp(—B/h)|1+O(h)]. The coefficient B was worked out in [38] 
as the on-shell Euclidean action of the bounce solution and the coefficient A was worked 
out in [39] to properly account for the quantum corrections t. In the very special case 
where the potential barrier is larger than the difference of energy density between false and 
true vacuums, a thin-wall approximation was proposed in [38] to evaluate the bounce ac- 
tion in a closed form consisting both contributions from the bubble interior and the bubble 
wall. The insight from thin-wall approximation provides us a physical picture of bubble 
nucleation that the released energy from spreading true vacuum into false vacuum goes to 
the acceleration of the bubble wall until a critical bubble with zero energy is formed. The 
size of such critical bubble was determined by the stationary point of the Euclidean action, 
which gives a rough cancelation between the energy from bubble interior and the energy 
from bubble wall. 

The seminal works |38, 39] were later extended to the case of vacuum decay in curved 
spacetime in [40] and thermal decay in flat spacetime in [41]. Although the presence of 
gravity does not change the general picture, it does lower the probability of vacuum decay 
by nucleating larger bubbles. In the extreme case of almost degenerated vacuums, gravity 
can even stabilize the false vacuum by preventing us from making a zero-energy bubble 
in the false vacuum, which in the absence of gravity can be simply done by expanding 
the bubble interior to balance the surface tension of bubble wall. When the temperature 
turns on, instead of an O(4)-symmetric bubble in Euclidean spacetime, an O(3)-symmetric 


'The precise form of A is not important in the study of GW from PTs because B dominates the evolution 
of decay probability though its exponential dependence. 


bubble solution that is periodic in time direction with period of inverse temperature T~! is 
expected, and the size of such critical bubble is determined by the maximum point of total 
energy. The tunneling rate was estimated in |[41, 42] as 


where the Euclidean action S3[(r), T], 


ogi Pe aay 
Ssle(r)] =40 | drf || —-) +V T|, (3.2) 
0 2 dr 
is estimated at the bounce profile ¢g(r) of equation-of-motion, 
2 
d'p 2dọ OV (3.3) 


dr? rdr ð 


The effective potential usually consists of tree-level potential, zero-temperature corrections 
and finite-temperature corrections including the daisy resummation. Recently, a new ther- 
mal resummation procedure was proposed in [43], which makes it possible to match theories 
to an EFT at finite temperature. Future works can be carried out along this direction. The 
naive strategy of solving above bounce equation is the shooting algorithm, which is illus- 
trated in Fig. 3. When multiple fields are involved at play, four approaches were introduced 
in the literatures, the early trial [44], the undamping/damping algorithm [45], the wildly 
adopted path deformation algorithm [|46], and the recently proposed semi-analytic pertur- 
bative approach [47]. 

The generalization of [38] in the case of non-zero cosmological constant allowed for 
both false and true vacuums was carried out in [48]. The vacuum decay is enhanced in 
the presence of gravity when the average of two vacuums is positive and suppressed when 
both vacuums are negative. The vacuum decay for a positive false vacuum but negative 
average of two vacuums is enhanced or suppressed when the gravitational effect becomes 
important or not. There is currently no satisfactory description of thermal decay in curved 
spacetime; however, we will clarify in the future work that the gravity correction only scales 
as k/(T R), and it is not important because the characteristic size of bubble R is much 
larger than the Planck length yk ~ 1//G for the characteristic scale of PT below the 
Planck scale. 


3.2 Bubble Expansion 


After the bubble is nucleated within false vacuum, the bubble wall will rapidly expand 
until approaching the speed of light [38] or colliding with other bubble walls. However, 


the realistic background is actually a thermal plasma full of relativistic particles 7, which 


In this review, we are only interested in the GWs from PTs that happened in radiation dominated 
era. For GWs from PTs that happened in inflation era and matter dominated era, please refer to [49-51] 
and [52], respectively. See also [53] for non-standard cosmology with an additional new component that 
redshifts faster than radiation. 
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Figure 3. The schematic illustration of shooting algorithm for the bounce equation, which is 
equivalent to a particle moving in the inverse potential with Hubble friction. The exit point of field 
profile exit presents the field value at the center of nucleated bubble, which can be found by the 
shooting algorithm with errors and trials of under/over shooting until an exact shooting is achieved. 
Here, the bubble profile is rescaled by the temperature Mg with respect to the bubble size rescaled 
by the mass scale Mp of effective potential. 


will interact with bubble walls by some form of frictions. The central problem of bubble 


expansion in plasma is to figure out the interconnections among following quantities: the 


—10- 


bubble wall velocity vw, the friction 7 on the wall from plasma, the strength a of PT that 
measure the released vacuum energy density with respect to the radiation energy density, 
and the efficiency factors kg and ky that measure the capability of transferring the liberated 
vacuum energy into the bubble wall expansion and the bulk fluid motion, respectively. These 
quantities serve as the interface between the theoretical models of particle physics and the 
simulations of bubble collisions. 

After the bubble is nucleated within thermal plasma, it starts stopping growing when 
the friction balances the net pressure from inside and outside of the bubble wall. If the 
pressure from outside of bubble wall is larger than the inside pressure while fluid velocity 
from outside of bubble wall is smaller than the inside velocity, the bubble wall behaves as 
deflagration, and the opposite definitions for detonation. Both deflagration and detonation 
are further classified as weak, Jouguet and strong types according to the fluid velocity 
from inside of bubble wall being smaller, equal and larger than the speed of sound. In the 
paper [54], it was proposed that with Jouguet condition, the velocity of bubble wall vw 
would be given by a simple formula [54] expressed in terms of the strength of PT a alone, 
so does the efficiency factor ky. The formula of Jouguet detonation was extensively used in 
the literatures at the early stage of studies of GWs from PTs despite the fact [55] that the 
Jouguet detonation can be unrealistic in the cosmological setup of PTs. Since the work [55], 
several parallel explorations have been considered as follows. 


e Beyond Chapman-Jouguet condition. Generalizing model-independent parametriza- 
tion equation [55] of friction, [56] was able to explore the full range of bubble wall 
velocity for both deflagrations and detonations, where the analytic approximations 
were found for both non-relativistic and ultra-relativistic wall velocities. In the state- 
of-art work of [57], it gives a unified picture and user-friendly fitting formulas for the 
dynamical regimes of bubble expansion. A different but more accurate approach was 
adopted in [58] right after [56] with microscopic considerations on the particle con- 
tents of specific models, which was revisited and improved in the subsequent work [59]. 
Future works can be carried out along this direction. 


e Criteria for runaway bubble walls. Apart from those stationary solutions of bubble 
wall with terminal velocity, there exists a runaway solution when the friction is too 
small to prevent the wall from approaching the speed of light. A simple criterion 
was found in [60] that the bubble walls will runaway if the effective potential in the 
true vacuum remains deeper than the false vacuum even after replacing the thermal 
potential by its second-order Taylor expansion term in the false vacuum, but not vice 
versa. This was reformulated in [57] simply by comparing the strength a of PTs with 
respect to a critical value a4. Hence, combined considerations of hydrodynamics and 
microphysics of specific models were later extended in [61] and [62] to the runaway 
regime. Future works can be carried out along this direction. 


e Reconciliations with baryogenesis. Large GWs signals from PTs require fast-moving 
bubble walls to collide with each other, while baryogenesis scenarios need slow-moving 
bubble wall to have an effective diffusion process. In [58], it also opens the possibility 
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3.3 


that can reconcile the baryogenesis with GWs from PTs in the presence of fermions 
with large Yukawa coupling and heavier stabilizing bosons. Later on, [63] pointed 
out that the relevant velocity responsible for the baryogenesis is the relative velocity 
between the wall and the plasma, which can be much smaller than the wall velocity 
when the strength of PTs becomes stronger. Therefore, it is possible to simultaneously 


generate large GWs from PTs without jeopardizing the effectiveness of baryogenesis. 


Bubble Percolation 


After bubble nucleation and bubble expansion, bubble percolation starts with colliding 
bubbles until PT is completed. When the initial size of nucleated bubble can be neglected, 
the duration of PT can be roughly estimated by the mean radius of bubbles at collisions, 


which is characterized by a single parameter 8. Both 6 and the previous mentioned a are 


evaluated at nucleation temperature T}, which is defined as the temperature at which the 


number of generated bubble per unit time per Hubble volume is of order one. The old 


picture of bubble percolation consists of two sources, namely, the colliding bubble walls and 


the turbulent motions of bulk fluids along with their associated magnetic field. However, 


the new picture of bubble percolation is added with an extra source from the sound waves of 


bulk motion as the result of bubble collision, which exist long after the bubble percolation. 


e Bubble collisions. It was first realized by Witten in [64] that the QCD PTs might 


leave us with detectable GWs from violent bubble collisions with their peak frequency 
characterized by the size of bubbles when they collided and with their peak amplitude 
estimated by the relative size of bubbles with respect to Hubble horizon at collision. 
Later on, this insight was generalized by Hogan in [65] to the case of EW PTs. Ina 
series of papers [66-69], the preliminary simulations were first implemented for bubble 
collisions to capture the general features of GWs from PTs. It is found in [66, 67] 
that remarkably the spectrum of GWs from simulating the collision of two vacuum 
bubbles in Minkowski space only depends upon the grossest features of the bubble 
collisions, namely, the strength a and duration 6 of PTs, which is similar to the result 
of simulation for hundreds of vacuum bubbles [68]. As we have shown in Fig. 4, the 
envelope approximation was proposed in [66-68] that the GWs are mainly generated 
from the uncollided envelopes of colliding bubble walls and any GWs from the overlap 
region can be neglected. The extension to the thermal bubble collisions was carried 
out later in the simulation [69], where the Jouguet detonation mode was explicitly 
used and the turbulent motion in fluid stirred by bubble collisions was appreciated 
in this case. The state-of-the-art results of the GWs spectrum from PTs due to 
colliding bubble walls were settled down in [70], which although disagreed with earlier 
study |71], they finally achieved consensus in [72]. Recently, under the thin-wall and 
envelope approximations in flat background, it was claimed in [73] that the GWs 
spectrum can be estimated analytically without need of simulations. Future works 
can be carried out along this direction from both theoretical and numerical points of 


view. 
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Figure 4. The schematic illustration of envelope approximation that the GWs are mainly generated 
from the uncollided envelopes of colliding bubble walls and any GWs from the overlap region can 
be neglected. 


e Turbulent MHD. The possibility of generating GWs from turbulent motion of bulk 
fluid was mentioned long time ago in [64] as remnant of bubble collisions, which 
was first estimated in [69] as Kolmogorov spectrum under quadrupole approximation. 
Apart from above turbulent velocity field, the fully ionized plasma could also gener- 
ate the turbulent magnetic field under turbulent motion, which itself is also a source 
of GWs. The early analysis |74, 75] involving MHD exhibits three problems [37]: 
large-scale problem (addressed in [76, 77]), time-evolution problem (persisted in [77] 
and addressed in |76]) and dispersion-relation problem (corrected in |76, 77|). The 
state-of-art result of the GWs spectrum from PTs due to non-helical MHD turbu- 
lence was analytically established in [78], where it ignored the circularly polarized 
GWs [79-85] from helical turbulence due to the macroscopic parity violation. The 
numerical simulations of relativistic MHD turbulence are certainly needed to make 
further confirmations for the above analytic results in the future. 


e Sound waves. The possibility of generating GWs from sound waves was originally 
pointed out in [65], although, it left forgotten for a long time until recent revelation 
in [86], where envelope approximation for modeling the colliding bubble walls is aban- 
doned and the GWs should be dominated by the overlapping sound waves in the bulk 
fluid. The breakthrough findings of [86] were further quantitatively understood in 


ae 


the updated simulations |87, 88] and theoretically modeled in [89], before which the 
possibility of having an inverse acoustic cascade was investigated in [90], suggesting 
the potentially very strong enhancement of the sound wave density at small wave 
number. Future works can be carried out along this direction for larger simulations 
with more bubbles. 


3.4 Gravitational Waves 


The GWs from the strong first-order PTs are guaranteed through above three processes: 
bubble nucleation, bubble expansion and bubble percolation. The bubble nucleation re- 
quires a potential barrier in order to tunnel from the false vacuum to the true vacuum. 
The bubble expansion requires fast enough moving bubble walls in order to have strong 
signals. The bubble percolation requires efficient collisions in order to dissipate released 
vacuum energy into bulk fluid motions. The fitting formulas of GWs spectrums from nu- 
merical simulations are well summarized in [36] and can be straightforwardly applied to 
the particle physics models, where the strength a and duration 8 evaluated at nucleation 
temperature T, along with the bubble wall velocity vw are obtained from the microscopic 
particle physics models, whereas the efficiency factors Kg and ky are approximated from the 
fitting formulas in [57]. There leaves one last discussion of which particle physics models 
can exhibit first-order PT so that GWs can be generated. 

The SM with phenomenological Higgs mechanism crosses over from the high-temperature 
phase to low-temperature phase [91] if the Higgs boson mass is larger than the W boson 
mass. To have GWs from first-order PTs, one has to go beyond SM (BSM) °. Here, we give 
an incomplete list of these models with explicit discussions on GWs, which will be revisited 
in the near future if we want to construct the reliable templates in order to extract the GW 


signals from the stochastic background. 


e Higher dimensional operators. The simplest example is to add a cubic term, which is 
usually expected in the supersymmetric extensions of SM in order to make a strong 
first-order PTs. By adopting the polynomial fitting formulae [93] for the bounce action 
from a general quartic potential with a cubic term, a semi-analytic calculation of the 
GW signal from the EWPT was carried out in [94]. The other important example 
is the dimension-six operator |95, 96], of which the GW signals from PTs were first 
analyzed in [97] and revisited in [98] with full scope of one-loop effective potential at 
finite temperature. In both examples, the bubble walls with detonations or run-away 
behavior have been considered in [99]. We produce in Fig. 5 the spectrum of GWs 
from PTs that originated from the dimension-six operator, which will be discussed in 
details in a future work [100]. 


e Additional scalar sectors. The simplest and extensively studied example is the gauge 
singlet scalar extension of SM [59, 101-107], which can be naturally fitted into previ- 
ously mentioned cubic term |94] and can be tested at future colliders by measuring the 


3Recently, it was argued in [92] that a thermally produced GWs signal could be generated in SM physics 
due to shear viscosity in the plasma. 
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Figure 5. The spectrum of GWs from PTs that originated from dimension-six operator with respect 
to the sensitivity curves of various configurations of eLISA project [37]. 


triple Higgs coupling precisely 4. The other important example is the charged scalar 
under SM gauge group, of which the simplest realization is the two-Higgs-doublet 
model (2HDM). The GW from 2HDM was preliminarily analyzed in [109] and further 
studied in [36] for the case of CP-conservation and recently revisited in [110] for the 
case of CP-violation. 


e Supersymmetric extensions. The capability of detecting GWs from PTs in Minimal Su- 
persymmetric Standard Model (MSSM) and Next-to-Minimal Supersymmetric Stan- 
dard Model (NMSSM) were first estimated in [111] and further explored in [112]. Both 
the PTs from MSSM [59] and NMSSM |97] were thought to be not strong enough to 
produce significant GW signals; however, the parameter space with strong first-order 
PTs has been identified for a modified version of MSSM |113] and a general version 
of NMSSM [62]. 


e Hidden dark sectors. The cosmological implications concerning with possible produc- 
tion of GWs from hidden dark sector were first explored in [114], and later discussed 
for light GeV scalar [115], vector thermal dark matter [116], UV-conformal dark sec- 
tor [117], SU(N) dark sectors with ny flavors [118], dark U(1) gauge complex scalar 
singlet [119], hidden sector with run-away bubble walls [120], and PTs involving suc- 
cessive hidden gauge symmetry breaking [121], dark matter asymmetry [122] and 


“Tt has been exemplified in [108] with dimension-six operator that the future space-based interferometers 
could probe the parameter spaces that are unreachable for the current particle colliders, but can have overlap 
with the future particle colliders, therefore a cross-check is possible. 
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two-step transition [123]. The GWs from first-order PTs could be a unique probe for 
these dark hidden sectors. 


e Other BSM extensions. The GWs from first order PTs were also analyzed in the case 
of extra dimensions |124-127], Peccei-Quinn PT [128], non-linear EWPT [129] and 
QCD PT [130-133]. 


4 The Gravitational Waves From Binary Systems 


GW sources can be divided into two categories, one is deterministic and the other is stochas- 
tic. For example, the primordial GW produced by the quantum fluctuations and by PTs 
during the early Universe belongs to the stochastic one. Due to the intrinsic random charac- 
ter, we cannot predict the waveform produced by the stochastic sources. The deterministic 
sources include two types. One type is predictable, while the other one is not. For example, 
the supernova is believed producing GWs. But the dynamics of supernova is too compli- 
cated to predict the related GW form. Another example is the foreground noise of LISA 
produced by white dwarf-white dwarf binaries in our milky way galaxy [134-136]. Due to 
the large number of such binaries, the combination of these GW signals makes the waveform 
unpredictable. In this section, we focus on the predictable GW sources. For such sources, 
we can construct theoretical model for gravitational waveform before GW detection exper- 
iment [137]. Then when the experiment data are ready, we can use matched filtering data 
analysis techniques to improve the GW detection sensitivity. And more, the theoretical 
model can be used to realize the astronomy detection of the GW sources through matched 
filtering [138, 139]. 

Binary systems are among the most important and major sources for GW detection 
projects including PTA, space-based laser interferometer detectors such as eLISA, Taiji 
and Tianqin projects, ground-based laser interferometer detectors such as Advanced LIGO, 
Advanced Virgo, Kagra and others. The GW frequency of binary system around merger 
state can be characterized by the ISCO (most Inner Stable Circular Orbit) frequency: 


1 


fisco = 63/297 M (4.1) 


where M is the total mass of the binary system, and we have used geometric unit with 
c= G= 1. When M = Mo, fisco ~ 2000Hz. At the same time, binary systems belong 
to predictable GW sources [139]. A typical example of the gravitational waveform and the 
related frequency are shown in Fig. 6. Note the time scale in the plot is different for inspiral 
stage and the merger/ringdown stage. The component of the binary system can be white 
dwarf, neutron star, black hole and even quark star. If gravitational theory beyond general 
relativity is true, other objects including axion star [140], gravastar [141] and other exotics 
might also be the component of binary system. 

Based on current understanding of stellar evolution, the final fate of a star may be 
white dwarf, neutron star or black hole which is determined by the mass of the star. Binary 
systems can be formed through capture process or many body interaction. Due to the 
gravitational slingshot effect involved in three-body interaction, binary systems can be 
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Figure 6. Typical example of the gravitational waveform and the corresponding frequency. In this 
example, the binary system is a binary black hole with equal mass. And the two black holes are 
spinless. The whole evolution process of a binary black hole system includes inspiral, merger and 
ringdown. But the boundary to distinguish these three stages is fuzzy. Compared to the inspiral 
stage, the merger and ringdown stages are much shorter. In this plot, we use different time scale 
for inspiral and merger /ringdown stage in order to make the merger/ringdown stage clearer. 


formed through many body interaction. Because a binary system is dynamically stable, 
binary systems are expected to be quite common in our Universe. 


The direct detection of GWs opens the GW astronomy. Through this new window to 
our Universe, GWs will not only bring us many new observations on kinds of objects in the 
Universe, but also give us many insights on fundamental physical theory. For example, black 
hole is a mystery from quantum gravity theory view. Especially, black hole presents us an 
information loss puzzle. Is it possible the black hole horizon be replaced with other objects 
such as firewall or anything else? Or any experiments present evidence that black hole 
horizon really exists? Unfortunately, one cannot get observational proof of the existence 
of a black hole event horizon based on EM waves |142]. In contrast, it is possible to use 
GW observation to give a direct proof of the existence of a black hole event horizon [143]. 
Currently, these promising aims can be most possibly achieved by predictable GW sources, 
especially binary systems. In order to realize these aims, we have to construct accurate 
enough gravitational waveform model for binary systems. There are three kinds of methods 
available to treat binary systems. They are PN approximation, numerical relativity and 
black hole perturbation method which are widely used for early inspiral stage, plunge and 
merger stage, and post merger stage, respectively. 
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4.1 Post-Newtonian Approximation 


Based on quadrupole formula [144], we can estimate the GWs related to a binary system: 
M 
h4 = —— 2° R*(1 + cos? 8) cos[2Q(t — r)], (4.2) 
M 02 p2 . 
hg = 74N R* cos 8 sin[2Q(t — r)], (4.3) 


where (r, 0, ġ) is the position of the observer with respect to the binary, M, R and Q are 
the total mass, the separation and the mutual orbiting frequency of the two components 
of the binary. Here h} and hx correspond to two polarization modes of GWs. These two 
modes are related to the dynamical metric form as 


ds? = —dt? + dr? + 77(1+h,)d6? + r° sin? 0(1 — h4 )dg? + 2r? sinOh,.dédd, (4.4) 


with (t,r,0,ġ) corresponding to the transverse-traceless gauge [145]. Here we need point 
out two points. First one is that the planer wave approximation has been taken in the 
above equation which is reasonable because the field point is very far away from the source. 
The second is that above estimation about h} and hx is only leading order of the PN 
approximation which is reasonable because the weak GW source moves slowly. For most 
realistic GW sources, these approximations will break down. Even for cases in which the 
approximation is reasonable, above approximation cannot satisfy the need of GW detec- 
tion [146, 147]. For example, when the two components of a binary system are some far 
away, they move slowly. So, this kind of situation corresponds to weak and moving slowly 
GW source. Although the above approximation is reasonable, the accuracy is not good 
enough. Then, more PN order corrections are needed to improve the accuracy [148]. 

The theoretical framework of the PN approximation for the processing of binary systems 
consists of two parts: the dynamics of binary components and the GWform. In order to 
construct the GW model, we need to solve the binary dynamics, and then put the solution 
into the waveform theory part to get the explicit GW model. In general, the dynamical 
equations of the PN approximation are a highly nonlinear system of ordinary differential 
equations. In that case, the numerical method has to be employed to get the solution. 
GW models in which the GWform is expressed in time domain include TaylorT2, TaylorT4 
and others [149]. Through stationary phase approximation, the post-Newtonian equations 
can be transformed into a frequency domain problem, and fortunately an approximate 
analytical expression can be obtained. Such a model is typically represented by the TaylorF2 
model. Corrected by the binary black hole spin dynamics, especially the effect of the orbital 
precession, the TaylorF2 model is modified and improved by single and double precession 
model [150, 151]. The TaylorT2 model is replaced by the X model when the elliptic orbit 
of the binary system is considered [152]. The theoretical model of the TaylorF2 model is 
extended to the elliptic orbit binary system including the post-circular model [153] and the 
enhanced post-circular model |154, 155]. 


4.2 Numerical Relativity 


Applying numerical methods to solve the Einstein equation is the topic of numerical rela- 
tivity. Currently, the numerical relativity can only deal with the problem of GW modeling 
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for the plunge and merger stage of binary systems. Because numerical relativity solves the 
Einstein equation without any approximation up to numerical error, this feature makes 
numerical relativity a universal tool to address kinds of GW sources. In Fig. 7, we show 
an example of the evolution process of a binary black hole system simulated by numeri- 
cal relativity. But the nature of diffeomorphism invariance of the general relativity brings 
numerical calculation special difficulties. The study of numerical relativity began in the 
60s of the last century [156], but the numerical instability made the code collapse after a 
few calculation steps. In 1990s, the construction of LIGO hardware strongly demanded the 
GW source modeling. So tens of universities and research institutions in USA and Europe 
jointly launched a Binary black hole Grand Challenge Project. However, the stability prob- 
lem of numerical relativity has not been solved by that project. Out of one’s expectation, 
the stability problem of numerical relativity was solved for the first time by Pretorius in 
2005 [157]. In the following 2006, Baker group and Campanelli group also independently 
solved the instability problem |158, 159]. Till now, there are more than ten numerical 
relativity groups around the world that have solved the stability problem, including the 
Princeton University, California Institute of Technology, University of Jena, Max Planck 


institute, the academy of mathematics and systems science, CAS [160] and others. 


It is worth pointing out that within kinds of stable numerical relativity methods, which 
tips and treatments are necessary and/or sufficient for stable computation are still an open 
question. From the viewpoint of the theory of partial differential equations, one can only 
apply the hyperbolicity analysis to linearized Einstein equation |161]. In practice, however, 
the success of Pretorius in 2005 depends much on the formalism of the Einstein equations 
and the adaptive mesh refinement technique. Regarding the formalism, the so-called BSSN 
equation and the generalized harmonic coordinate equation have been widely used in numer- 
ical relativity. Adaptive mesh refinement is a very effective method to deal with multiscale 
problems. At present, the parallel adaptive mesh refinement codes developed especially for 
Einstein equations include the BAM code (developed by Bruegmann), AMSS-NCKU code 
(developed by Zhoujian Cao) [162], PAMR code (developed by Pretorius) and the Carpet 
code (developed by Schnetter). 


Besides the stability problem, the major issues involved in numerical relativity are 
accuracy and efficiency. As for the formalism factor, more and more investigations show 
that Z4c and CCZ4 formalism is better than the BSSN formalism [163-165]. As an example 
of partial differential equation formalism of Einstein equations, we show Z4c formalism as 
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Figure 7. Typical evolution of binary black hole simulated by numerical relativity. The yellow 
region corresponds to the highly curved region which is roughly the position of black holes. The 
four snapshots correspond to inspiral, merger and ringdown stages. This result is simulated by 
AMSS-NCKU software. 


following, which is proposed first time in 3D form in [164]: 
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The unknown functions which need to be solved numerically are listed in the left-hand side 
of the above equations. For the meaning of the notations used in the above equations, 
we refer our readers to [164]. Spectral method, finite difference method and finite element 
method are three categories of numerical methods for partial differential equations. The vast 
majority of numerical relativity groups take the finite differential method. The SpEC code of 
California Institute of Technology uses spectral method. On the other hand, the application 
of finite element method into numerical relativity is still very few [|166]. The spectral 
method has good computational efficiency due to its exponential convergence. However, the 
characteristics of its global data exchange limits its ability of parallel computing. The finite 
difference method working with domain decomposition algorithm can achieve good parallel 
computing efficiency. In order to deal with the multiscale characteristics of astrophysics, 
the refinement of the multilayer data structure is essential in the finite difference method. 
However, this method limits the number of cores for parallel computation by a single data 
layer, thus limits the strong parallel scalability. In contrast, the finite element method 
can combine the exponential convergence of the spectral method and the high parallel 
scalability of the difference method. Therefore, in principle, it is expected that the finite 
element calculation of the Einstein equations can achieve good strong parallel scalability. 
However, it is still an issue to construct the weak form of the Einstein equations, and to 
realize large-scale scientific computation [137, 139]. 

Now the major challenges to numerical relativity placed by GW source modeling are 
the calculation of binary systems with mass ratio more than 100 to 1 and the well-converged 
simulation of binary systems including neutron stars [138]. 


4.3 Black Hole Perturbation 


After the merger of the binary components, the system will ringdown and finally settle 
down to a Kerr black hole. If the two components are white dwarf and/or neutron star, 
the final remnant could be a stable neutron star. Here, we only concern with the case 
with Kerr black hole as the final product. Then the ringdown stage can be described 
by perturbation of a Kerr black hole. The black hole perturbation theory was pioneered 
by Regge, Wheeler and Zerilli [167, 168] in which perturbation around a Schwarzschild 
black hole was considered. Teukolsky investigated the perturbation of a Kerr black hole 
on spacetime curvature level [169, 170]. The former is called the Regge-Wheeler-Zerilli 
equation, and the latter is called the Teukolsky equation. The perturbation method used 
by Teukolsky is also valid to Schwarzschild black hole. For Schwarzschild black hole, Regge- 
Wheeler-Zerilli method and Teukolsky method is equivalent [171]. In the early 70s of the 
last century, some scholars have begun to use the Regge-Wheeler-Zerilli equation to study 
the associated GW for a test particle falling into a black hole [172-174]. Because the 
Teukolsky equation admits spin of a black hole, it is valid more extensively than the Regge- 
Wheeler-Zerilli equation. In recent years, there have been many works using the Teukolsky 
equation to calculate the GWs. 

The Teukolsky equation can be solved by the PN approximation method (note that 
the Teukolsky equation is the target in question here, instead of Einstein equation itself 
like what described in the above subsection) or through numerical method. The PN ap- 
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proximation method was mainly developed by Mano, Suzuki, Takasugi and others. These 
authors developed some hypergeometric functions and Coulomb wave functions to approx- 
imate the homogeneous solution of the Teukolsky equation [175, 176]. Later on, Fujita and 
his coworkers applied PN approximation method to these analytical solutions to calculate 
the GWform and the related energy flux for extreme mass ratio binary systems. Currently, 
22PN accuracy for Schwarzschild black hole and 11PN accuracy for Kerr black hole have 
been achieved |177, 178]. However, such accuracy still do not yet satisfy the requirement 
of GW detection experiment [179]. 

The numerical methods to solve the Teukolsky equation can be divided into two cate- 
gories: frequency domain method and time domain method. The frequency domain method 
divides the original Teukolsky equation into two ordinary differential equations through vari- 
able separating method and the Fourier expansion [169]. Regarding the radial equation, 
we can transform it to Sasaki-Nakamura equation before solving it numerically [180]. The 
Teukolsky equation can only provide the energy flux and the GWs. In order to provide 
the orbit of the test particle, Hughes and his coworkers applied adiabatic approximation 
method to the geodesic equation to treat the inspiral [181]. Now this kind of approxima- 
tion has become an important method to treat binary systems with extreme mass ratio. 
Different from the usual finite difference numerical method, Fujita and Tagoshi used the 
hypergeometric function and the Coulomb wave function to expand the original equations 
which gives a quicker and more accurate numerical method [182]. The numerical method 
in time domain needs to solve a set of 2+1 partial differential equations. At present, one 
mainly uses the finite difference method to solve the problem (but see [183] for the finite 
element method). Due to the extreme mass ratio of the involved binary system, the related 
GW signal may last several years. Although the computation required by Teukolsky equa- 
tion is much simpler than numerical relativity, how to improve the computational efficiency 
of the Teukolsky equation is also a great challenge for GW source modeling. 


4.4 Cosmological Probes 


In 1986, Schutz showed that from the GWs one can infer the Hubble constant by using the 
fact that the GWs from the binary systems encode the absolute distance information |184]. 
The inspiraling and merging compact binaries consisting of neutron stars and black holes 
can be considered as standard candles, or “standard sirens”. The name of “siren” is due 
to the fact that the GW detectors are omni-directional and detect coherently the phase 
of the wave, which makes them in many ways more like microphones for sound than like 
conventional telescopes. For an expanding Universe, the chirp waveform of the GW can be 
generalized to the cosmological case by multiplying all masses by the factor 1 + z and the 
physical distance D can be replaced by the luminosity distance dz |185, 186]. In fact, the 
waveform h produced by compact binary inspirals during the inspiral phase is theoretically 
well described by the analytical solution: 


a (MY (AL) cosesinl (4.11) 
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which is valid at the lowest (Newtonian) order for the “cross” GW polarization (the same 
expression holds for the “plus” polarization with a difference dependence on the orientation 
of the binary’s orbital plane v). Here, Me(z) is the (redshifted) chirp mass, f the GW 
frequency at the observer and ®(f) its phase. For the standard sirens, the luminosity 
distance dz(z) is the most important parameter entering the waveform (4.11). Parameter 
estimation over the observed GW signal (4.11) can directly yield the value of the luminosity 
distance of the GW source, together with an uncertainty due to the detector noise. Once the 
signal is detected and characterized, the luminosity distance of the source can be extracted. 
This indicates that one can measure the luminosity distance directly without the need of 
the cosmic distance ladder: standard sirens are “self-calibrating”. However, it is not possible 
to infer the redshift z from the GW signal of the binary of black holes because all of the 
observed parameters such as the masses and distance are redshifted by the same factor 
1+ z. To use the distance measurement for cosmography, one has to obtain the redshifts 
of the GW sources by some independent methods. 

Before considering how to obtain the redshift information, one should think about how 
accurately the distance can be measured. The performance of the standard siren is limited 
by several effects. Firstly, the intrinsic measurement uncertainty in the amplitude of the 
detector’s response is simply the inverse of the signal-to-noise ratio (SNR) [187, 188], which 
is related to the detector sensitivity. Second, the weak gravitational lensing will distort 
the measurement of the luminosity distance, producing a magnification or demagnification 
on the order of a few percent [189-192]. Furthermore, the largest contribution of the 
uncertainty comes from the limited direction and source orientation sensitivity, and there 
is a large correlation between the distance, the sky position and orientation of the source. 
Thus, a network of detectors is needed to measure the position and orientation of the 
binary and to break the degeneracy among these parameters. By observing a simultaneous 
detection of a beamed EM signal, one can determine the sky position accurately and also 
will help to improve the measurements of the distance. The EM signal, called the EM 
counterpart, is also an important issue to the GW on cosmography. 

The most traditional way to obtain the redshift of a GW event is through an accom- 
panying EM signal, the EM counterpart. The binary merger of an NS with either an NS 
(BNS) or BH (BHNS) is hypothesized to be the progenitor of a short and intense burst of 
y-rays, a so-called SGRB [193]. An EM counterpart like SGRB can provide the redshift 
information if the host galaxy of the event can be pinpointed. Moreover, SGRBs are likely 
to be strongly beamed phenomena [194], which allow one to constrain the inclination of the 
compact binary system, breaking the distance-inclination degeneracy. The GWs with short 
y-rays bursts or other EM counterparts as the standard sirens have been studied in various 
papers (see [187, 188, 190, 191, 195-210], and references therein). For example, Nissanke 
et.al. [202] used MCMC method and found that a network of advanced LIGO detectors can 
constrain the Hubble constant to a 5% accuracy. In Ref. [203], the authors demonstrate 
that with 1000 GW events detected by ground-based Einstein Telescope ° it is possible to 
constrain the Hubble constant and dark matter energy parameter up to Ahg ~ 5 x 1073 
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and AQ, ~ 0.02 using Fisher matrix approach, which was also found by Cai et al. [188] 
using MCMC method. Furthermore, Cai et al. used the Gaussian Process method and 
found that the equation of state of the dark energy can be constrained to Aw(z) ~ 0.03 in 
the low redshift region, which gives a better constraint than Ref. [203]. For the space-based 
detector LISA °, the expansion of the Universe and interacting dark energy have also been 
studied by Refs. |208, 210]. 

There also exist other methods to infer the redshift information of the GW events, such 
as galaxy catalog proposed by Schutz [184], neutron star mass distribution [195, 211], and 
the tidal deformation of neutron stars |212, 213]. Measuring the redshift associated with 
a GW event is one of the biggest challenges in the future, see [214]; however, for a new 
GW standard sirens probe without redshift information by utilizing those BH binaries to 
be a tracer of the large-scale structure [215]. In addition, the spin of BH can also help 
us estimate the GW’s parameters. GWs as the standard sirens to probe the cosmological 
parameters can provide an independent and complementary alternative to current experi- 
ments. It is expected that combining GW data with other astronomical observations such 
as supernovae, and adopting a better data analysis approach, the cosmological parameters 


could be constrained more precisely than the current situation. 


5 Conclusions 


The direct detection of GWs by LIGO initiates a new era of GW astronomy and GW 
cosmology. The GW physics is not only related to gravitational physics, but also closely 
related to particle physics, cosmology and astrophysics. The GWs provide us a new powerful 
tool to reveal various secrets of the nature. Indeed, a lot of relevant papers have appeared 
since the announcement of the direct detection of GWs. In this paper, we have briefly 
introduced three kinds of GW sources and relevant physics. They are GWs produced 
during inflation and preheating in the early Universe, from cosmic PTs and dynamics of 
compact binary systems, respectively. We also have discussed in a simple way the GWs as 
standard siren in the evolution of the Universe. Due to the limitation of space, we are not 
able to discuss all aspects of GW physics, but only focus on some main issues. Of course, 
it is also impossible to list a complete list of the references, quite probably some important 


references are missed here, we should apologize for this. 
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